function [out]=varDecHor(GG,RR,ZZ,SDX,NVDec)
% Called from momentsDraws.m
% out.varDecHorObs: [NZ NX NVDEC], Variance Decomposition Across Horizons for Observables 
% out.varDecHorStates: [NS NX NVDEC]  Decomposition of the *whole* state
% vector across horizons 
[ny nx]=size(RR);
nz=size(ZZ,1); 
Gnew=eye(ny);
out.varDecHorObs=zeros(nz,nx,NVDec);
out.varDecHorStates=zeros(ny,nx,NVDec);
vzero=RR*(SDX')*SDX*(RR');
vcum=zeros(nz,nz,NVDec,nx);
vcum_st=zeros(ny,ny,NVDec,nx);
for kk=1:NVDec;
    % ====================
    % Whole Variance
    % ====================
    if kk==1;
        vcfull_st=vzero;
        vcfull=ZZ*vcfull_st*(ZZ');
    else
        Gnew=GG*Gnew ;
        vcfull_st=Gnew*vzero*(Gnew')+vcfull_st;
        vcfull=ZZ*Gnew*vzero*(Gnew')*(ZZ')+vcfull;
    end;
    deno_temp=diag(vcfull);
    deno_st=diag(vcfull_st);
    if kk==1;
        for jj=1:nx
            et=( SDX(jj,:)' )*SDX(jj,:);
            vcum_st(:,:,1,jj)=(RR*et*(RR'));
            vcum(:,:,1,jj)=ZZ*((RR*et*(RR')))*(ZZ');
            out.varDecHorObs(:,jj,kk)=(diag(vcum(:,:,1,jj)))./deno_temp;
            out.varDecHorStates(:,jj,kk)=(diag(vcum_st(:,:,1,jj)))./deno_st;
        end;
    else
        for jj=1:nx;
            et=( SDX(jj,:)' )*SDX(jj,:);
            vcum(:,:,kk,jj)=ZZ*(Gnew*RR*et*(RR')*(Gnew'))*(ZZ')+vcum(:,:,kk-1,jj);
            vcum_st(:,:,kk,jj)=(Gnew*RR*et*(RR')*(Gnew'))+vcum_st(:,:,kk-1,jj);
            out.varDecHorObs(:,jj,kk)=(diag(vcum(:,:,kk,jj)))./deno_temp;
            out.varDecHorStates(:,jj,kk)=(diag(vcum_st(:,:,kk,jj)))./deno_st;
        end
    end
    if any( sum( squeeze(out.varDecHorObs(:,:,kk)) , 2 ) < 0.99 )
        error('Variance Decompositions Observables are innacurate');
    end
    
    if any( sum( squeeze(out.varDecHorStates(:,:,kk)) , 2 ) < 0.99 )
        error('Variance Decompositions Observables are innacurate');
    end
end

end

